In [412]:
def timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x, simulation=0,title=''):  
    plt.figure(figsize=(10, 6))
    plt.subplot(311)
    plt.plot(x, C[:,simulation], alpha = .5)
    plt.ylabel('Coral Node Count')

    plt.subplot(312)
    plt.ylabel('CN Index')
    plt.xlabel('Time')
    plt.plot(x, CN[:,simulation], alpha = .5)
    
    plt.subplot(313)
    plt.ylabel('CNr Index')
    plt.xlabel('Time')
    plt.plot(x, CNr[:,simulation], alpha = .5)

def initialFinal9(types=types, C=C, CN=CN, CNr=CNr, simulation=1, 
                 timesteps=[0,1,2,3,4,5,6,7,8,9,10,11,12]):
    plt.figure(figsize=(15, 12))

    #colors = ['pink', 'lightgreen','darkgreen']
    #levels = [0, 1, 2]
    #cmap, norm = clt.from_levels_and_colors(levels=levels, colors=colors, extend='max')
    
    plt.subplot(331).set_title('-C-  ' + str(C[timesteps[0],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[0],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[0],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[0],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    plt.subplot(332).set_title('-C-  ' + str(C[timesteps[1],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[1],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[1],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[1],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    plt.subplot(333).set_title('-C-  ' + str(C[timesteps[2],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[2],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[2],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[2],:,simulation],10), 
               cmap=cmap, norm=norm)  
    
    plt.subplot(334).set_title('-C-  ' + str(C[timesteps[3],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[3],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[0],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[3],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    plt.subplot(335).set_title('-C-  ' + str(C[timesteps[4],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[4],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[4],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[4],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    
    plt.subplot(336).set_title('-C-  ' + str(C[timesteps[5],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[5],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[5],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[5],:,simulation],10), 
               cmap=cmap, norm=norm) 
    
    
    plt.subplot(337).set_title('-C-  ' + str(C[timesteps[6],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[6],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[6],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[6],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    
    plt.subplot(338).set_title('-C-  ' + str(C[timesteps[7],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[7],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[7],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[7],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    
    plt.subplot(339).set_title('-C-  ' + str(C[timesteps[8],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[8],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[8],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[8],:,simulation],10), 
               cmap=cmap, norm=norm)
    plt.show()

def initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=1, 
                 timesteps=[0,1,2]):
    plt.figure(figsize=(15, 4))

    #colors = ['pink', 'lightgreen','darkgreen']
    #levels = [0, 1, 2]
    #cmap, norm = clt.from_levels_and_colors(levels=levels, colors=colors, extend='max')
    
    plt.subplot(131).set_title('-C-  ' + str(C[timesteps[0],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[0],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[0],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[0],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    plt.subplot(132).set_title('-C-  ' + str(C[timesteps[1],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[1],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[1],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[1],:,simulation],10), 
               cmap=cmap, norm=norm)
    
    plt.subplot(133).set_title('-C-  ' + str(C[timesteps[2],simulation])+
                              '  -CN-  ' + str(round(CN[timesteps[2],simulation],2))+
                              '  -CNr-  ' + str(round(CNr[timesteps[2],simulation],2)),
                              fontweight="bold")
    plt.imshow(shaper(types[timesteps[2],:,simulation],10), 
               cmap=cmap, norm=norm)  
    plt.show()

Welcome to the Coral Model Testing and Exploration Notebook!

In this notebook, we take you through examples of running our model the various way to plot the model's outputs.

We begin with loading the necessary libraries

In [1]:
from coralModelTracking import Organism, Reef, Ocean

import tools as tl
import numpy as np

import matplotlib.pyplot as plt
import matplotlib.colors as clt

from mpl_toolkits.axes_grid1.inset_locator import InsetPosition

Next, we set out model parameters

In [2]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100

## Grid shapes
rows, columns = 10, 10

## Inital Coral/Algae Percentages
coralPercent = .5
algaePercent = .25

## Parameters
r=1.0
d=.4 #.2
a=.2
g=.9
y=.75 #.9

threshold = 1.45 #extent of neighborhood recognition

NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent

Now we create and run multiple reefs

In [3]:
Pacific = Ocean()

for s in range(0,NumberOfSimulations):

    Moorea = Reef()
    count = 0
    for i in range(0,rows):        
        for j in range(0,columns):
            U = np.random.choice([0,1,2],
                                 p=[coralPercent, turfPercent, algaePercent])
            node = Organism(type=U, location=[i,j], ID=count)
            Moorea.append(node)
            count = count + 1
            
    Moorea.generateGraph(threshold)

    for n in range(0,NumberOfTimesteps):
        Moorea.roll(r=r, d=d, a=a, g=g, y=y, dt=dt) 
        
    Pacific.append(Moorea)

The Products we now can play with

  • Pacific.simulation[1].coralNodeCount
  • Pacific.simulation[1].coralNeighborCount

Functions to pull information from these objects

Can pull keys and values via:

  • x = list(dictionary.keys()), and
  • y = list(dictionary.values())

In running the coralModel, we are interested in capturing the CN index: the average number of coral neighbors for each coral node.

Single Reef Example, looking at Simulation 1

In [4]:
example = Pacific.simulation[1].coralNodeCount
x = list(example.keys())            ## the time steps
coralCount = list(example.values()) ## the number of coral nodes

plt.figure(figsize=(15, 5))
plt.plot(x,coralCount)
Out[4]:
[<matplotlib.lines.Line2D at 0x1a18285240>]
In [5]:
example = Pacific.simulation[1].coralNeighborCount
x = list(example.keys())
CN = np.array(list(example.values()))/np.array(coralCount) ## taking the average of the neighborhood count

plt.figure(figsize=(15, 5))
plt.plot(x,CN)
plt.plot(x,np.repeat(4, NumberOfTimesteps))
Out[5]:
[<matplotlib.lines.Line2D at 0x1a18435c18>]

Each node can have up to 8 neighbors. If the average is above 4 (the orange line above), this means on average, each coral node is surrounded by more coral than all other benthic types combined.

In [6]:
example = Pacific.simulation[1].coralNeighborCount
x = coralCount
y = CN

plt.figure(figsize=(15, 5))
plt.scatter(x,y)
Out[6]:
<matplotlib.collections.PathCollection at 0x1a185d0470>

Looking at the Average across the various simulations

In [7]:
def dictToNumpy(dictionary):
    output = list(dictionary.values())
    output = np.array(output)
    return(output)

multsimCoralCounts = np.array([dictToNumpy(Pacific.simulation[i].coralNodeCount) 
                               for i,val in enumerate(Pacific.simulation)]).transpose()

multsimNeighborCounts = np.array([dictToNumpy(Pacific.simulation[i].coralNeighborCount) 
                               for i,val in enumerate(Pacific.simulation)]).transpose()

Coral Node Count: 7 simulations and average of 100 simulations

In [10]:
x = list(Pacific.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts[:,1:7], alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[10]:
[<matplotlib.lines.Line2D at 0x1a18746828>]

CN index: 7 simulations and average of 100 simulation

In [11]:
plt.figure(figsize=(15, 5))
CN = multsimNeighborCounts[:,1:7]/multsimCoralCounts[:,1:7]
CNmean = (multsimNeighborCounts/multsimCoralCounts).mean(axis=1)
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[11]:
[<matplotlib.lines.Line2D at 0x102252828>]

Simulation with Maximized Coral Neighborhood Blob

In [12]:
## distance to center mask
center = (round(rows/2),round(columns/2))
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center) 
                 for i in range(0,rows) 
                 for j in range(0,columns)]).reshape(rows,columns)
plt.figure(figsize=(15, 5))
plt.imshow(distanceGrid)
Out[12]:
<matplotlib.image.AxesImage at 0x1a18e9c6d8>

Checking Initial Reef

In [64]:
coralNodeLocations = (np.where(distanceGrid < np.median(distanceGrid)))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n]) 
                      for n in range(0,len(coralNodeLocations[0]))]
### Sudo Making of Initial Reef
grid = np.zeros((rows,columns))
for i in range(0,rows):        
    for j in range(0,columns):
        U = np.random.choice([1,2],
                             p=[.5, .5])
        if (i,j) in coralNodeLocations:
            U = 0
        grid[i,j] = U
plt.figure(figsize=(15, 5))
plt.imshow(grid)
Out[64]:
<matplotlib.image.AxesImage at 0x1a218f28d0>

Model set-up/run with initial coral blob condition, coral at 50%

In [15]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100

## Grid shapes
rows, columns = 10, 10

## Inital Coral/Algae Percentages
coralPercent = .5
algaePercent = .25

## Parameters
r=1.0
d=.4 #.2
a=.2
g=.9
y=.75 #.9

threshold = 1.45 #extent of neighborhood recognition

NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent

## Specific Coral Nodes
center = (round(rows/2),round(columns/2))
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center) 
                 for i in range(0,rows) 
                 for j in range(0,columns)]).reshape(rows,columns)
coralNodeLocations = (np.where(distanceGrid < np.median(distanceGrid)))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n]) 
                      for n in range(0,len(coralNodeLocations[0]))]
In [16]:
## Running Model
Resilient = Ocean()
for s in range(0,NumberOfSimulations):

    optimalReef = Reef()
    count = 0
    for i in range(0,rows):        
        for j in range(0,columns):
            U = np.random.choice([1,2],
                                 p=[.5, .5])
            if (i,j) in coralNodeLocations:
                U = 0
            node = Organism(type=U, location=[i,j], ID=count)
            optimalReef.append(node)
            count = count + 1            
    optimalReef.generateGraph(threshold)

    for n in range(0,NumberOfTimesteps):
        optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt) 
        
    Resilient.append(optimalReef)

Coral Node Count: 7 simulations and average of 100 simulations

In [65]:
multsimCoralCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNodeCount) 
                               for i,val in enumerate(Pacific.simulation)]).transpose()

multsimNeighborCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNeighborCount) 
                               for i,val in enumerate(Pacific.simulation)]).transpose()
##plotting
x = list(Resilient.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts[:,1:7], alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[65]:
[<matplotlib.lines.Line2D at 0x1a21902eb8>]

CN index: 7 simulations and average of 100 simulations

In [19]:
plt.figure(figsize=(15, 5))
CN = multsimNeighborCounts[:,1:7]/multsimCoralCounts[:,1:7]
CNmean = (multsimNeighborCounts/multsimCoralCounts).mean(axis=1)
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[19]:
[<matplotlib.lines.Line2D at 0x1a1a66da58>]

Model set-up/run with initial coral blob condition, coral at 33%

In [20]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100

## Grid shapes
rows, columns = 10, 10

## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33

## Parameters
r=1.0
d=.4 #.2
a=.2
g=.9
y=.75 #.9

threshold = 1.45 #extent of neighborhood recognition

NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent


## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center) 
                 for i in range(0,rows) 
                 for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]

coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n]) 
                      for n in range(0,len(coralNodeLocations[0]))]
In [21]:
%%time
Resilient = Ocean()

for s in range(0,NumberOfSimulations):

    optimalReef = Reef()
    count = 0
    for i in range(0,rows):        
        for j in range(0,columns):
            U = np.random.choice([1,2],
                                 p=[.5, .5])
            if (i,j) in coralNodeLocations:
                U = 0
            node = Organism(type=U, location=[i,j], ID=count)
            optimalReef.append(node)
            count = count + 1            
    optimalReef.generateGraph(threshold)

    for n in range(0,NumberOfTimesteps):
        optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt) 
        
    Resilient.append(optimalReef)
CPU times: user 1min 10s, sys: 590 ms, total: 1min 11s
Wall time: 1min 11s

Coral Node Count: 7 simulations and average of 100 simulations

In [23]:
multsimCoralCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNodeCount) 
                               for i,val in enumerate(Pacific.simulation)]).transpose()

multsimNeighborCounts = np.array([dictToNumpy(Resilient.simulation[i].coralNeighborCount) 
                               for i,val in enumerate(Pacific.simulation)]).transpose()
x = list(Resilient.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts[:,1:7], alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[23]:
[<matplotlib.lines.Line2D at 0x1a190e8710>]

CN index: 7 simulations and average of 100 simulations

In [24]:
plt.figure(figsize=(15, 5))
CN = multsimNeighborCounts[:,1:7]/multsimCoralCounts[:,1:7]
CNmean = (multsimNeighborCounts/multsimCoralCounts).mean(axis=1)
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[24]:
[<matplotlib.lines.Line2D at 0x1a1a500198>]

G parameter = .5, Model set-up/run with initial coral blob condition , coral at 33%

In [80]:
## Simulations and Number of runs
NumberOfSimulations = 100
dt=.1
tf=100

## Grid shapes
rows, columns = 10, 10

## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33

## Parameters
r=1.0
d=.4 #.2
a=.2
g=.5
y=.75 #.9

threshold = 1.45 #extent of neighborhood recognition

NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent


## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center) 
                 for i in range(0,rows) 
                 for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]

coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n]) 
                      for n in range(0,len(coralNodeLocations[0]))]
In [81]:
%%time
LowGrazing = Ocean()

for s in range(0,NumberOfSimulations):

    optimalReef = Reef()
    count = 0
    for i in range(0,rows):        
        for j in range(0,columns):
            U = np.random.choice([1,2],
                                 p=[.5, .5])
            if (i,j) in coralNodeLocations:
                U = 0
            node = Organism(type=U, location=[i,j], ID=count)
            optimalReef.append(node)
            count = count + 1            
    optimalReef.generateGraph(threshold)

    for n in range(0,NumberOfTimesteps):
        optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt) 
        
    LowGrazing.append(optimalReef)
CPU times: user 1min 1s, sys: 175 ms, total: 1min 1s
Wall time: 1min 2s

Coral Node Count: 100 simulations and average of 100 simulations

In [82]:
multsimCoralCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNodeCount) 
                               for i,val in enumerate(LowGrazing.simulation)]).transpose()

multsimNeighborCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNeighborCount) 
                               for i,val in enumerate(LowGrazing.simulation)]).transpose()
x = list(LowGrazing.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts, alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[82]:
[<matplotlib.lines.Line2D at 0x1a24c91748>]

CN index: 100 simulations and average of 100 simulations

In [83]:
a = multsimNeighborCounts
b = multsimCoralCounts

CN = np.divide(a, b, out=np.zeros_like(a), where=b!=0)
CNmean = CN.mean(axis=1)

plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[83]:
[<matplotlib.lines.Line2D at 0x1a24fb4d68>]

Same as above but with more time steps, testing with less simulations

In [84]:
## Simulations and Number of runs
NumberOfSimulations = 10
dt=.1
tf=1000

## Grid shapes
rows, columns = 10, 10

## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33

## Parameters
r=1.0
d=.4 #.2
a=.2
g=.5
y=.75 #.9

threshold = 1.45 #extent of neighborhood recognition

NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent


## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center) 
                 for i in range(0,rows) 
                 for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]

coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n]) 
                      for n in range(0,len(coralNodeLocations[0]))]
In [85]:
%%time
LowGrazing = Ocean()

for s in range(0,NumberOfSimulations):

    optimalReef = Reef()
    count = 0
    for i in range(0,rows):        
        for j in range(0,columns):
            U = np.random.choice([1,2],
                                 p=[.5, .5])
            if (i,j) in coralNodeLocations:
                U = 0
            node = Organism(type=U, location=[i,j], ID=count)
            optimalReef.append(node)
            count = count + 1            
    optimalReef.generateGraph(threshold)

    for n in range(0,NumberOfTimesteps):
        optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt) 
        
    LowGrazing.append(optimalReef)
CPU times: user 1min 2s, sys: 363 ms, total: 1min 2s
Wall time: 1min 3s
In [86]:
multsimCoralCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNodeCount) 
                               for i,val in enumerate(LowGrazing.simulation)]).transpose()

multsimNeighborCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNeighborCount) 
                               for i,val in enumerate(LowGrazing.simulation)]).transpose()

Coral Node Count: 10 simulations of 10,000 timesteps, and the average of the 10 simulations

In [87]:
x = list(LowGrazing.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts, alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[87]:
[<matplotlib.lines.Line2D at 0x1a253119b0>]

^^ For this one above it would be interesting to explore whether the spatial configuration for the two top lines allowed for coral success.

CN index: 10 simulations of 10,000 time steps, and the average of the 10 simulations

In [88]:
a = multsimNeighborCounts
b = multsimCoralCounts

CN = np.divide(a, b, out=np.zeros_like(a), where=b!=0)
CNmean = CN.mean(axis=1)

plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[88]:
[<matplotlib.lines.Line2D at 0x1a21315e80>]

G parameter = .2, more time steps, testing with less simulations

In [89]:
## Simulations and Number of runs
NumberOfSimulations = 10
dt=.1
tf=1000

## Grid shapes
rows, columns = 10, 10

## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33

## Parameters
r=1.0
d=.4 #.2
a=.2
g=.2
y=.75 #.9

threshold = 1.45 #extent of neighborhood recognition

NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent


## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center) 
                 for i in range(0,rows) 
                 for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]

coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n]) 
                      for n in range(0,len(coralNodeLocations[0]))]
In [90]:
%%time
LowGrazing = Ocean()

for s in range(0,NumberOfSimulations):

    optimalReef = Reef()
    count = 0
    for i in range(0,rows):        
        for j in range(0,columns):
            U = np.random.choice([1,2],
                                 p=[.5, .5])
            if (i,j) in coralNodeLocations:
                U = 0
            node = Organism(type=U, location=[i,j], ID=count)
            optimalReef.append(node)
            count = count + 1            
    optimalReef.generateGraph(threshold)

    for n in range(0,NumberOfTimesteps):
        optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt) 
        
    LowGrazing.append(optimalReef)
CPU times: user 57.9 s, sys: 369 ms, total: 58.2 s
Wall time: 58.7 s
In [91]:
multsimCoralCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNodeCount) 
                               for i,val in enumerate(LowGrazing.simulation)]).transpose()

multsimNeighborCounts = np.array([dictToNumpy(LowGrazing.simulation[i].coralNeighborCount) 
                               for i,val in enumerate(LowGrazing.simulation)]).transpose()

Coral Node Count: 10 simulations of 10,000 timesteps, and the average of the 10 simulations. Very grazing (g=.2)

In [92]:
x = list(LowGrazing.simulation[1].coralNodeCount.keys())
plt.figure(figsize=(15, 5))
plt.plot(x,multsimCoralCounts, alpha=0.4)
plt.plot(x,multsimCoralCounts.mean(axis=1), 'black')
Out[92]:
[<matplotlib.lines.Line2D at 0x1a254c8048>]

Coral Node Count: 10 simulations of 10,000 timesteps, and the average of the 10 simulations. Very low grazing (g=.2)

In [93]:
a = multsimNeighborCounts
b = multsimCoralCounts

CN = np.divide(a, b, out=np.zeros_like(a), where=b!=0)
CNmean = CN.mean(axis=1)

plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[93]:
[<matplotlib.lines.Line2D at 0x1a1b155898>]

-- 30 SIMULATIONS, GRAZING .5, SPATIAL PLOTS ---

In [266]:
## Simulations and Number of runs
NumberOfSimulations = 30
dt=.1
tf=1000

## Grid shapes
rows, columns = 10, 10

## Inital Coral/Algae Percentages
coralPercent = .33
algaePercent = .33

## Parameters
r=1.0
d=.4 #.2
a=.2
g=.5
y=.75 #.9

threshold = 1.45 #extent of neighborhood recognition

NumberOfTimesteps = int(tf/dt)
NumberOfNodes = rows * columns
turfPercent = 1 - coralPercent - algaePercent

## Storing Spatial Grid, may slow down run
#types = np.zeros((NumberOfTimesteps, NumberOfNodes, NumberOfSimulations))

## Distance Prep
distanceGrid = np.array([Reef.distance([i+.5,j+.5], center) 
                 for i in range(0,rows) 
                 for j in range(0,columns)]).reshape(rows,columns)
maxCoralDist = np.concatenate(distanceGrid, axis=None)
maxCoralDist = np.sort(maxCoralDist)[33]

coralNodeLocations = (np.where(distanceGrid < maxCoralDist))
coralNodeLocations = [(coralNodeLocations[0][n],coralNodeLocations[1][n]) 
                      for n in range(0,len(coralNodeLocations[0]))]
In [267]:
%%time
#Switching = Ocean()

for s in range(0,NumberOfSimulations):

    optimalReef = Reef()
    count = 0
    for i in range(0,rows):        
        for j in range(0,columns):
            U = np.random.choice([1,2],
                                 p=[.5, .5])
            if (i,j) in coralNodeLocations:
                U = 0
            node = Organism(type=U, location=[i,j], ID=count)
            optimalReef.append(node)
            count = count + 1            
    optimalReef.generateGraph(threshold)

    for n in range(0,NumberOfTimesteps):
        
        for i,val in enumerate(optimalReef.nodes):
            types[n,i,s] = optimalReef.nodes[i].type
        optimalReef.roll(r=r, d=d, a=a, g=g, y=y, dt=dt) 
        
    Switching.append(optimalReef)
CPU times: user 3min 26s, sys: 764 ms, total: 3min 27s
Wall time: 3min 28s
In [274]:
# Save types, i.e. all the grid progression data
#for i in range(0,len(Switching.simulation)):
#    np.savetxt("coralModelTrackingOutputs /modelOutput_switching"+ str(i)+".csv", 
#               np.reshape(types[:,:,i], (-1, rows)), delimiter=",")
In [441]:
#import pickle
#path = 'coralModelTrackingOutputs/'
#output_name = 'grazing5'
#outfile = open(path+output_name, "wb")       #open pickle jar
#pickle.dump(Switching, outfile)           #put contents into jar
#outfile.close()                         #close the jar

#PacificReturns = pickle.loads(open(path+output_name, "rb").read())

Coral Count over time, 30 Simulations

In [376]:
x = list(Switching.simulation[1].coralNodeCount.keys())

multsimCoralCounts = np.array([dictToNumpy(Switching.simulation[i].coralNodeCount) 
                               for i,val in enumerate(Switching.simulation)]).transpose()

multsimNeighborCounts = np.array([dictToNumpy(Switching.simulation[i].coralNeighborCount) 
                               for i,val in enumerate(Switching.simulation)]).transpose()
C = multsimCoralCounts
N = multsimNeighborCounts

plt.figure(figsize=(15, 5))
plt.plot(x,C, alpha=0.4)
plt.plot(x,C.mean(axis=1), 'black')
Out[376]:
[<matplotlib.lines.Line2D at 0x3fb7c121d0>]

CN index over time, 30 Simulations

In [351]:
CN = np.divide(N, C, out=np.zeros_like(N), where=C!=0)
CNmean = CN.mean(axis=1)

plt.figure(figsize=(15, 5))
plt.plot(x, CN, alpha=0.4)
plt.plot(x, CNmean, 'black')
Out[351]:
[<matplotlib.lines.Line2D at 0x3f88482f60>]

CNr index over time, 30 Simulations

In [352]:
CNr = np.divide(CN, C, out=np.zeros_like(a), where=b!=0)
CNrmean = CNr.mean(axis=1)

plt.figure(figsize=(15, 5))
plt.plot(x, CNr, alpha=0.4)
plt.plot(x, CNrmean, 'black')
Out[352]:
[<matplotlib.lines.Line2D at 0x3f89447ba8>]

Conclusions from exploring each of the simulations:

  • Lives On (simulation:time at which relative steadyness is reached)
  • simulation 7:400, 14:1100, 25:900

  • Mid Crash (simulation:crashing time)

  • simulation 2:2300-3000, 22:3000, 28:3000

  • EarlyMid Crash (simulation:crashing time)

  • simulation 8:1800, 13:1700, 17:1800, 21, 23

  • Ambiguous early crashers (simulation:crashing time)

  • simulation 0:1000, 1:300, 3:400, 4:300, 5:800, 6, 9, 10, 11, 12, 15, 16, 18, 19, 20, 24, 26, 27, 29

One task for tomorrow: run these grid configurations again and see if how their outputs vary (maybe this is where to include the variance)?

Plots where coral prevails - simulations 7, 14, 25

C, CN, CNr Time series example simulation 14

In [444]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=14, 
                  title="Makes it for a while, and then crashes. grazing =.5")

Spatial evolutions of simulations 7, 14, and 25

In [464]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=7, 
                 timesteps=[0,100,500])
In [465]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=14, 
                 timesteps=[0,1000,2000])
In [466]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=25, 
                 timesteps=[0,700,1000])

Plots from mid-crashes, simulations 2, 22, 28

C, CN, CNr Time series example simulation 8

In [417]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=2, 
                  title="Makes it for a while, and then crashes. grazing =.5")

Spatial evolutions of simulations 2, 22, and 28

In [461]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=2, 
                 timesteps=[0,2000,3000])
In [462]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=22, 
                 timesteps=[0,2000,3000])
In [463]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=28, 
                 timesteps=[0,2000,3000])

Plots from Early-Mid Crashes, simulations 8, 13, 17

C, CN, CNr time series example simulation 8

In [428]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=8, 
                  title="Makes it for a while, and then crashes. grazing =.5")

Spatial evolutions of simulations 8, 13, 17

In [467]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=8, 
                 timesteps=[0,1000,2000])
In [468]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=13, 
                 timesteps=[0,1000,2000])
In [469]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=17, 
                 timesteps=[0,1000,2000])

Plots from Early Crashes, simulations 6, 9, 10

C, CN, CNr time series example simulation 6

In [432]:
timeSeriesPlotter(C=C,CN=CN,CNr=CNr,x=x,simulation=6, 
                  title="Makes it for a while, and then crashes. grazing =.5")

Spatial evolutions of simulations 6, 9, 10

In [470]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=6, 
                 timesteps=[0,100,300])
In [471]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=9, 
                 timesteps=[0,100,300])
In [472]:
initialFinal3(types=types, C=C, CN=CN, CNr=CNr, simulation=10, 
                 timesteps=[0,100,300])

Looking at plots with same C but heading towards different states

Time Series simulation 14

In [487]:
timeSeriesPlotter(C=C[0:3000,:],CN=CN[0:3000,:],
                  CNr=CNr[0:3000,:],x=x[0:3000],simulation=14)
In [459]:
initialFinal9(types=types, C=C, CN=CN, CNr=CNr, simulation=14, 
                 timesteps=[0,10,75,
                           90,110,125,
                           230,240,250])

Time Series simulation 2 (crashes)

In [488]:
timeSeriesPlotter(C=C[0:3000,:],CN=CN[0:3000,:],
                  CNr=CNr[0:3000,:],x=x[0:3000],simulation=2)
In [460]:
initialFinal9(types=types, C=C, CN=CN, CNr=CNr, simulation=2, 
                 timesteps=[0,10,75,
                           90,110,125,
                           230,240,250])

Exploration of weird CNr vs CN effect

Coral Dominating Simulations

In [490]:
simulation = 7
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[490]:
<matplotlib.collections.PathCollection at 0x3f72d68240>
In [491]:
simulation = 14
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[491]:
<matplotlib.collections.PathCollection at 0x3fa8587978>
In [492]:
simulation = 25
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[492]:
<matplotlib.collections.PathCollection at 0x3f6f1389e8>

Mid Crashes

In [493]:
simulation = 2
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[493]:
<matplotlib.collections.PathCollection at 0x3fa9037208>
In [495]:
simulation = 22
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[495]:
<matplotlib.collections.PathCollection at 0x3f746aae80>
In [496]:
simulation = 28
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[496]:
<matplotlib.collections.PathCollection at 0x3f691e6198>

Coral Early-Mid Crashes

In [497]:
simulation = 8
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[497]:
<matplotlib.collections.PathCollection at 0x2cc7cb21d0>
In [498]:
simulation = 13
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[498]:
<matplotlib.collections.PathCollection at 0x3f69693710>
In [499]:
simulation = 17
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[499]:
<matplotlib.collections.PathCollection at 0x3faa6e5208>

Coral Early Crashes

In [501]:
simulation = 6
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[501]:
<matplotlib.collections.PathCollection at 0x3f7382aba8>
In [502]:
simulation = 9
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[502]:
<matplotlib.collections.PathCollection at 0x3faa821c50>
In [503]:
simulation = 10
plt.figure(figsize=(15, 5))
plt.scatter(CN[:,simulation], CNr[:,simulation], alpha=0.1)
Out[503]:
<matplotlib.collections.PathCollection at 0x3f6fc517b8>
In [ ]: